function [VelDironeVelDataMask VelDironeVelDataAverage VelDironeVelDataUncertainty...
    flowRateVsTime averageUncertainty averageVelocity avgVelocityVsTime peakVelocityVsTime CrossSectionArea FlowRate]...
    = FlowStats(VelDirOneVelData,ARGUS_ROI, BW, PixelSpacingRows, PixelSpacingColumns, t, ROInumber)

% Find stats -------------------------------------------------------------%
% calculate uncertainty at each pixel ------------------------------------%
[Rsize Csize Phase] = size(VelDirOneVelData);
VelDironeVelDataUncertainty = zeros(Rsize,Csize,Phase);
VelDironeVelDataMask        = zeros(Rsize,Csize,Phase);
VelDironeVelDataAverage     = zeros(Rsize,Csize,Phase);

h = waitbar(0,['Calculating Velocities for ROI' num2str(ROInumber,1)]);

if (ARGUS_ROI == 1) | (ARGUS_ROI ==3) | (ARGUS_ROI ==4)
    tempTimeIndex = Phase;
    CrossSectionArea(1:tempTimeIndex) = 0;

else
    tempTimeIndex = 1; % don't recurse BW if not temporal
    CrossSectionArea = 0;
end

for Rindex = 1:1:Rsize
    waitbar(Rindex/Rsize)
    
    for Cindex = 1:1:Csize
        for TimeIndex = 1:1:tempTimeIndex
           
            if BW(Rindex,Cindex,TimeIndex)
                VelDironeVelDataMask(Rindex,Cindex,TimeIndex) =...
                    VelDirOneVelData(Rindex,Cindex,TimeIndex);
                VelDironeVelDataAverage(Rindex,Cindex)     = ...
                    mean(VelDirOneVelData(Rindex,Cindex,:));
                VelDironeVelDataUncertainty(Rindex,Cindex) = ...
                    std(VelDirOneVelData(Rindex,Cindex,:)) / ...
                    sqrt(length(VelDirOneVelData(Rindex,Cindex,:)));
                CrossSectionArea(TimeIndex) = ...
                    CrossSectionArea(TimeIndex) + PixelSpacingRows * ...
                    PixelSpacingColumns;
            end
        end
    end
end
close(h);

h = waitbar(0,['Calculating Averages for ROI' num2str(ROInumber,1)]);
% Averages . . . 
dumIndex = 1;
for C = 1:1:Csize
    waitbar(C/Csize)
       
    for R = 1:1:Rsize
        if BW(R,C,1)  % need to fix this for varying ROIs
            VelDirOneAverageArray( dumIndex ) = ...
                VelDironeVelDataAverage(R,C);
            VelDirOneUncertaintyArray( dumIndex ) = ...
                VelDironeVelDataUncertainty(R,C);
            dumIndex = dumIndex + 1;
        end
    end
end

close(h)

peakVelocity    = max(VelDirOneAverageArray);
averageVelocity = mean(VelDirOneAverageArray);

peakUncertainty    = max(VelDirOneUncertaintyArray);
averageUncertainty = mean(VelDirOneUncertaintyArray);

%Calculate the temporal mean and max velocity versus time
for TimeIndex = 1:t
    peakVelocityVsTime(TimeIndex)    = max(max(VelDironeVelDataMask(:,:,TimeIndex)));
    minVelocityVsTime(TimeIndex)     = min(min(VelDironeVelDataMask(:,:,TimeIndex)));

    if abs(minVelocityVsTime(TimeIndex)) > peakVelocityVsTime(TimeIndex)
        peakVelocityVsTime(TimeIndex) = minVelocityVsTime(TimeIndex);
    end

    avgVelocityVsTime(TimeIndex) = mean(mean(nonzeros(VelDironeVelDataMask(:,:,TimeIndex))));
    flowRateVsTime(TimeIndex) = avgVelocityVsTime(TimeIndex)*CrossSectionArea(TimeIndex)*.0006;

end

% Average Flow Rate
FlowRate = averageVelocity * (mean(CrossSectionArea/100)); %mL/s
FlowRate = FlowRate * .06; %L/min